function GSSOut = gs_gridgene(GSS)
%gs_gridgene - XY-grid generation
%
%       GSSOut = gs_gridgene(GSS)
%
% This function generates a regular grid on which the strain field can be 
% computed.
% These new fields are added to GSS leading to GSSOut 
%       GridReport
%       XGridStep
%       YGridStep
%       XGridNumber
%       YGridNumber
%       DownLeftCorner
%       TopRightCorner
%
% Default: XGridStep (dx) and YGridStep (dy) are the 50% of standard 
% deviations of x and y coordinate dispersion respectively, with the 
% constrain that the not lower among them is an integer multiple of the 
% other (i.e., if dx <= dy, it is dy = m*dx for a positive integer m). 

% A. Pesci & G. Teza, 2004, 2005, 2022.

format long e;

e = GSS.East;
n = GSS.North;
npe = GSS.IncludedPoints;
IMS = GSS.RefImageStruct;

xmax = max(e); 
xmin = min(e);
ymax = max(n); 
ymin = min(n);

%--- Points included/not included in the calculations ---
ei = e;
ni = n;
ei(~npe) = NaN;
ni(~npe) = NaN;
en = e;
nn = n;
en(npe) = NaN;
nn(npe) = NaN;

% --- grid preliminary individuation ---
xx = sort(e);
yy = sort(n);
dx = std(diff(xx));
dy = std(diff(yy));
%-------------------------------------------------------------
[dx,dy,~] = multy(dx,dy);  % Ancillary function: see below
%-------------------------------------------------------------
grx = xmin-3*dx:dx:xmax+3*dx;
gry = ymin-3*dy:dy:ymax+3*dy;
nx  = length(grx);
ny  = length(gry);
fvx = grx(1);  
fvy = gry(1);
lvx = grx(nx); 
lvy = gry(ny);  

GSSOut = GSS;
GSSOut.XGridStep = dx;
GSSOut.YGridStep = dy;
GSSOut.XGridNumber = nx;
GSSOut.YGridNumber = ny;
GSSOut.DownLeftCorner = [fvx fvy];
GSSOut.TopRightCorner = [lvx lvy];

nconty = 0;
while nconty == 0
    
    [X,Y] = meshgrid(grx,gry);

    grep = gs_gridReport(GSSOut);
    
    % visualization -------------------------------
    figk = figure;
    set(figk,'Units','normalized','OuterPosition',[0 0 1 1]);
    plot(ei,ni,'xk',en,nn,'xm');
    legend('Point included','Point not used','Location','NorthWestOutside');
    hold on
    plot(X,Y,'g','HandleVisibility','off');
    plot(X',Y','g','HandleVisibility','off');   % in order to draw the grid
    xlabel('EAST'); 
    ylabel('NORTH'); 
    title('GRID AND DATA POINTS');
    ALA = axis;
    gs_ReadBGImage(IMS,ALA);
    axis equal
    hold off
    %--------------------------------------------------
    
    lb = uicontrol(figk,'Style','listbox',...
        'string', {...
            'ACCEPT THE GRID AS IT IS',...
            'CHANGE GRID STEPS ONLY',...
            'CHANGE GRID LIMITS ONLY',...
            'CHANGE GRID STEPS AND LIMITS',...
            'NO GRID GENERATION'},...
        'FontSize',14,...
        'Units','normalized','Position', [0.01 0.05 0.3 0.15],...
        'Callback', 'uiresume(gcbf)');
    txt = uicontrol(figk,'Style','text',...
        'Units','normalized','Position',[0.01 0.4 0.3 0.3],...
        'String',{grep},...
        'FontSize',14);
    uiwait(figk);
    lbv = lb.Value;
    delete(lb);
    delete(txt);
    
    if lbv == 5
        GSSOut = GSS;
        return
    end
    
    if (lbv == 2) || (lbv == 4)
        %-------------------------------------------------------
        [dx,dy] = new_d(dx,dy);
        %---------- Ancillary functions: see below -------------
		[dx,dy,~] = multy(dx,dy);  
        %-------------------------------------------------------
        fvxp = fvx; % grid step change -> vertices change!
        fvyp = fvy;  
        lvxp = lvx; 
        lvyp = lvy;
    end
    if (lbv == 3)||(lbv == 4)
         lbgl = uicontrol(figk,'Style','listbox',...
            'string',...
                {sprintf('NEW GRID LIMITS FROM FIGURE %d (DOWN-LEFT, TOP-RIGHT)',figk.Number),...
                    'NEW GRID LIMITS FROM KEYBOARD',...
                    'SKIP GRID LIMITS CHANGE'},...
            'FontSize',14,...
            'Units','normalized','Position',[0.05 0.05 0.40 0.15],...
            'Callback', 'uiresume(gcbf)');
         uiwait(figk);
         lbglv = lbgl.Value;
         delete(lbgl);
         if lbglv == 1
             %--------------------------------------------------------- 
             [fvxp,fvyp] = new_l(fvx,fvy,figk,'DOWN-LEFT'); 
             %----------- Ancillary function: see below ---------------
             [lvxp,lvyp] = new_l(lvx,lvy,figk,'TOP-RIGHT',fvxp,fvyp);
             %---------------------------------------------------------
         elseif lbglv == 2
             %---------------------------------------------------------
             [fvxp,fvyp,lvxp,lvyp] = new_lik(fvx,fvy,lvx,lvy);
             %---------------------------------------------------------
         else
             lbv = 1;
         end
    end
    if lbv > 1
        grx = fvxp:dx:lvxp;
		gry = fvyp:dy:lvyp;
		nx  = length(grx);
		ny  = length(gry);
		fvx = grx(1);  
        fvy = gry(1);
        lvx = grx(nx); 
        lvy = gry(ny);
        
    else
        nconty = 1;
    end
    
    grep = gs_gridReport(GSSOut);
    GSSOut.GridReport = grep;
    GSSOut.XGridStep = dx;
    GSSOut.YGridStep = dy;
    GSSOut.XGridNumber = nx;
    GSSOut.YGridNumber = ny;
    GSSOut.DownLeftCorner = [fvx fvy];
    GSSOut.TopRightCorner = [lvx lvy];
    
    close(figk);
    
end


%% ANCILLARY FUNCTIONS

function [dx,dy,nm] = multy(dx,dy)
% This function makes the not lower number between dx and dy an integer 
% multiple of the other. 
% The integer nm is the ratio.
% If min(dx,dy) is very near to zero, dx = dy = max(dx,dy) is used. 
dmin = min(dx,dy);
dmax = max(dx,dy);
if abs(dmin) < 1000*eps
    nm = 1;
    dx = dmax;
    dy = dmax;
else
    nm = ceil((dmax-dmin+1000*eps)/dmin);
    % (ceil(x) is the nearest integer to x towards plus infinity) 
    if abs(dmin-dx) < 1000*eps
        dy = nm*dx;
    else
        dx = nm*dy;
    end
end

function rd = nu_dd(r,nd)
% For the representation of a number as '%nip.ndf', in which nip is the 
% number of digits of its integer part and nf is the number of decimals.
srd = ['%' num2str(length(num2str(fix(r)))) '.' num2str(nd) 'f'];
     % (The result is the format for the representation)
rd = sprintf(sprintf('%s',srd),r);

function [dxn,dyn] = new_d(dx,dy)
% Interactive choice of new grid steps.
prompt = {'Grid step along x-axis','Grid step along y-axis'};
dxs = num2str(dx);
dys = num2str(dy);
defans = {dxs,dys};
fields = {'dx','dy'};
acquis = inputdlg(prompt,'GRID STEPS MODIFICATION',1,defans,'on');
if ~isempty(acquis)              %see if user hit cancel
    acq = cell2struct(acquis,fields);
    dxn = str2double(acq.dx);      %convert string to number
    dyn = str2double(acq.dy);
else
    dxn = dx;
    dyn = dy;
end
if isempty(dxn)
    dxn = dx; 
end
if isempty(dyn)
    dyn = dy; 
end

function [vx,vy] = new_l(vx,vy,figk,namevert,vxmi,vymi)
% Interactive choice of vertex 'namevert', with coordinates (vx,vy), picked 
% on figure k and with the condition that vx > vxmi, vy > vymi.
% If vxmi and vxma are unassigned, the default values vxmi = -Inf, vymi = -Inf 
% are used.
if nargin < 6
    vxmi = -Inf; 
    vymi = -Inf; 
end
figure(figk);
nocont = 0;
while nocont == 0
    [vx,vy] = ginput(1);
    if (vx > vxmi) && (vy >= vymi)
        strvx = nu_dd(vx,2); 
        strvy = nu_dd(vy,2);
        [indx,tf] = listdlg(...
            'PromptString',...
                {sprintf('THE NEW SELECTED %s CORNER IS (%s,%s)',...
                    namevert,strvx,strvy)'},...
            'SelectionMode','single',...
            'ListString',{'ACCEPT','REJECT'});
        if tf
            menok = indx;
        else
            menok = 2;
        end
        if menok == 1
            nocont = 1; 
        end
    end
end

function [fvxn,fvyn,lvxn,lvyn] = new_lik(fvx,fvy,lvx,lvy)
% Interactive choice of the vertices using the keyboard
defansy = {num2str(fvx), num2str(fvy),...
    num2str(lvx), num2str(lvy)};
fields = {'fvx','fvy','lvx','lvy'};
nocont = 0;
nit = 1;
while nocont == 0
    if nit == 1
        prompty = {'Down-left corner, east coordinate',...
            'Down-left corner, north coordinate',...
            'Top-right corner, east coordinate',...
            'Top-right corner, north coordinate'};
    else
        prompty = {sprintf('%s\n\n%s',...
            'ERROR: THE CORNERS ARE DOWN-LEFT AND TOP-RIGHT',...
            'Down-left corner, east coordinate'),...
            'Down-left corner, north coordinate',...
            'Top-right corner, east coordinate',...
            'Top-right corner, north coordinate'};
    end
    acquis = inputdlg(prompty,'GRID LIMITS',1,defansy);
    if ~isempty(acquis)
        acq  = cell2struct(acquis,fields);
        fvxn = str2double(acq.fvx);
        fvyn = str2double(acq.fvy);
        lvxn = str2double(acq.lvx);
        lvyn = str2double(acq.lvy);
        if (fvxn < lvxn) && (fvyn < lvyn)
            nocont = 1;
        end
    else
        fvxn = fvx; 
        fvyn = fvy;
        lvxn = lvx; 
        lvyn = lvy;
        nocont = 1;
    end
    nit = nit+1;
end